Relativistic Turbulence: A Long Way from Preheating to Equilibrium 
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We study, both numerically and analytically, the development of equilibrium after preheating. We 
show that the process is characterised by the appearance of Kolmogorov spectra and the evolution 
towards thermal equilibrium follows self-similar dynamics. Simplified kinetic theory gives values for 
all characteristic exponents which are close to what is observed in lattice simulations. The resulting 
time for thermalization is long, and temperature at thermalization is low, T ~ 100 eV in the simple 
A$ 4 inflationary model. Our results allow a straightforward generalization to realistic models. 



Introduction. The dynamics of equilibration and 
thermalization of field theories is of interest for various 
reasons. In high-energy physics understanding of these 
processes is crucial for applications to heavy ion collisions 
and to reheating of the early universe after inflation. In- 
flation solves the flatness and the horizon problems of 
the standard big bang cosmology and provides a calcu- 
lable mechanism by which initial density perturbations 
were generated [1]. At the end of inflation the Universe 
was in a vacuum-like state. In the process of decay of 
this state and subsequent thermalization (reheating) the 
matter content of the universe is created. It was realized 
recently that the initial stage of reheating, dubbed pre- 
heating [2], is a fast, explosive process. This initial stage 
by now is well understood [3-7]. Strong and fast ampli- 
fication of fluctuation fields at low momenta may lead 
to various interesting physical effects, like non-thermal 
phase transitions [8], peculiar baryogenesis [9], genera- 
tion of high-frequency gravitational waves [10], etc. 

Understanding of the subsequent stages of reheating 
and thermalization processes and calculation of the final 
equilibrium temperature is important for various appli- 
cations, most notably baryogenesis and the problem of 
over-abundant gravitino production in supergravity mod- 
els [11]. Thermalization of field theories was discussed 
already, see e.g. Refs. [12]. However, at present the pro- 
cess of thermalization after preheating is still far away 
from being well understood and developed. The problem 
is that at the preheating stage the occupation numbers 
are very large, of order of the inverse coupling constant. 
In addition, in many models the zero mode does not de- 
cay completely. Therefore, a simple kinetic approach is 
not applicable. 

Fortunately, the description in terms of classical field 
theory is valid in this situation [3], and the process of 
preheating, as well as subsequent thermalization, can be 
studied on a lattice. In this paper we adopt this ap- 



proach. Our goal is to integrate the system on a lattice 
sufficiently accurately and sufficiently far in time to be 
able to see generic features, and possibly to the stage, at 
which the kinetic description becomes a good approxima- 
tion scheme. Lattice studies of thermalization, similar to 
ours, were done in Ref. [13]. Several generic rules of ther- 
malization were formulated, like the early equipartition of 
energy between coupled fields. However, the problem is 
very complicated and there are other unanswered impor- 
tant questions like what is the final thermalization tem- 
perature, at what stage the kinetic description becomes 
valid, what is the functional form of particle distributions 
during the thermalization stage, etc. 

For our study we use a higher accuracy, improved ver- 
sion of the LATTICEEASY code [14]. We show that the 
distribution functions follow a self-similar evolution re- 
lated to the turbulent transport of wave energy. This 
property enables us to estimate the physical reheating 
temperature, which turns out to be very low. The con- 
cept should be rather model independent since typical 
ranges of particle momenta at preheating and in ther- 
mal equilibrium are widely separated. However, in this 
letter we will restrict our numerical integration (but not 
the discussion) to the "minimal" inflationary model, the 
massless A<E> 4 -theory. 

The Model. With conformal coupling to gravity and 
after a rescaling of the field, cp = $a, where a(t) is the 
cosmological scale factor, the equation of motion in co- 
moving coordinates describes a (/^-theory in Minkowski 
space-time, 

n<p + A^ 3 = 0. (1) 

At the end of inflation the field is homogeneous, ip = 
tpo(t). Later on fluctuations develop, but the homoge- 
neous component of the field, which corresponds to the 
zero momentum in the Fourier decomposition, may be 
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dynamically important and is referred to as the "zero- 
mode." In such situations it is convenient to make a 
further rescaling of the field, cj) = <p/<Po(t>o)i an d of the 
space-time coordinates, — > y/Xipo(t )x tJ - , which trans- 
forms the equation of motion (1) into dimensionless and 
parameter free form, 



U(f> + 4? = o . 



(2) 



Here to corresponds to the initial moment of time (end 
of inflation), and in what follows we denote dimension- 
less time as r. With this rescaling the initial condition 
for the zero- mode oscillations is 0o( T o) = 1- All model 
dependence on the coupling constant A and on the ini- 
tial amplitude of the field oscillations now is encoded in 
the initial conditions for the small (vacuum) fluctuations 
of the field with non-zero momenta [3]. The physical 
normalization of the inflationary model corresponds to a 
dimensionful initial amplitude of y>o(*o) ~ 0.3Mpi and a 
coupling constant A ~ 10 -13 [1]. The re-parametrization 
property of the system allows to chose a larger value of 
A for numerical simulations. We have used A = 10~ 8 . 

Numerical Procedure and Results. We use a 3-D cubic 
lattice with periodic boundary conditions. The finite- 
differences scheme that was used is 2nd order in time and 
4-th order in space. The results displayed here are taken 
from a simulation with 256 3 sites and a physical box size 
L = 14%. With this box size the infrared modes which 
belong to the resonance band are still well represented, 
while the ultraviolet lattice cut-off is sufficiently far away 
from the occupied modes, such that the particle spectra 
are not distorted even at late times. We have studied the 
dependence of our results on the lattice- and the box size 
to avoid lattice artifacts. Various quantities are measured 
and monitored both in configuration space (zero mode, 
4>a = (0), and the variance, var(<fi) = (4> 2 ) — and 
in the Fourier space. Using fourier transformed fields 
we first define the wave amplitudes (which correspond to 
annihilation operators in the quantum problem), 



a(k) 



LOk<Pj: + i(f>k 
(2tt) 3 / 2 V2^ 

The effective frequency ujk = 



(3) 



+ rn e ff is determined 

by the effective mass rn\^ = 3A(0 2 ) and the inverse fc 2 ) 
of the lattice Laplacian. In our numerical scheme the 
latter is given by 



kl = b-* 



5 8 I 
^ . g - 3 cos(fcfci) + - cos(26fci) ) . (4) 

ie{i,2,3} '* ' 



Here b = 2ir/L — 1/7 is the lattice constant. Making 
use of afc, we calculate various correlators, n(k) = (a'a), 
a(k) = (aa), (aWetet), etc. The first one, which corre- 
sponds to the particle occupation numbers, is of prime 
interest. 
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FIG. 1. Amplitude of the zero-mode oscillations, cj> , and 
variance of the field fluctuations as functions of time r. 
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FIG. 2. Occupation numbers as function of k<f> Q 
t = 100, 400, 2500, 5000, 10000. 
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We begin the discussion of our numerical results with 
the evolution of the zero-mode and the variance of the 
field, which are shown in Fig. 1. Initially we see an 
exponentially fast transfer of the zero-mode energy into 
fluctuations during preheating (up to t~ 300). It is fol- 
lowed by a long and slow relaxation phase. In this regime 
(t > 1500) the amplitude of the zero mode oscillations 
decreases according to ~ r -1 / 3 , the variance of the field 
(averaged over high-frequency oscillations) drops accord- 
ing to ~ t~ 2 / 5 . This is consistent with previous results 
[3] . In addition we find that in this regime the zero-mode 
is in a non-trivial dynamical equilibrium with the bath 
of highly occupied modes: when the zero- mode is artifi- 
cially removed, it is recreated on a short time-scale (Bosc 
condensation). 

At early times the distribution functions of particles 
over momenta, see Fig. 2, have peaky structure. The 
first peak which corresponds to the parametric resonance 
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is initially at the theoretically predicted value of A: ~ 1.27 
[5]. Later (r > 1500) the spectra become smooth and at 
small k approach a power-law, rik ~ k~ s , where s fluctu- 
ates in the range of 1.5 — 1.7, depending on time and the 
range of k where it is fitted. This power law clearly dif- 
fers from the the classical thermal equilibrium, rift ~ ■ 
It is followed by the exponential cut-off, whose position 
monotenously shifts with time towards higher k. Pump- 
ing of energy from the zero mode stays effective all the 
times (note a small bump in the particle distributions in 
Fig. 2 at k ~ 1). It corresponds to the annihilation of 
four condensed particles into two quanta. Rescattering 
of two particles into two particles is also effective. One 
of the two can belong to the zero-mode condensate either 
in initial or in the final two particle state. We also can 
see in Fig. 2 that in the power-law region rik is a func- 
tion of k/4> only, where </> ( T ) represents the amplitude 
of the zero-mode at time r. This effect can be related 
to the above described dynamical equilibrium between 
zero-mode and the bath of particles. Indeed, going from 
Eq. (1) to Eq. (2) we can rescale by the current ampli- 
tude of the zero-mode. 

The picture presented in Fig. 2 at late times resembles 
stationary Kolmogorov turbulence. It appears as such 
due to rescaling of momenta by the amplitude of the zero- 
mode, but in fact in the present model the turbulence can 
not be stationary because the amplitude of the zero-mode 
(i.e. the strength of the source of turbulence) decreases. 
Further examination of Fig 2 suggests that the evolution 
of particle spectra may be self-similar. We have tried 
therefore the following anzatz 



n(k,T) = T- q n Q (kT- p ). 



(5) 



Spectra rescaled at several moments of time by the rela- 
tion inverse to Eq. (5) are shown in Fig. 3. We have found 
that the evolution is indeed self-similar with q sa 3.5p and 
p » 1/5. 

Discussion. Here we discuss the question whether a 
simple kinetic theory gives predictions for turbulence and 
self-similarity exponents in agreement with the lattice 
calculations. Our lattice study of higher order correla- 
tors, like (a'a'aa), shows that the field distribution is 
very close to Gaussian, see also [12,13]. This facilitates 
the use of the kinetic approach. On the other hand we 
have found that the magnitude of a(k) is of order of a 
few percent compared to n(k), and it is even larger in 
the region of resonant momenta. This means that the 
strict kinetic approach should include a(k). Nevertheless 
we neglect these effects and write a kinetic equation in 
a simple form rik — Ik, where the collision integral for a 
TO-particle interaction is given by 



dn k U k F\n] 



(6) 



In d spatial dimensions the integration measure d£lk is 
given by m — 1 integrations over d-dimensional Fourier 
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FIG. 3. On the right hand side we plot the wave energy per 
decade found in lattice integration. On the left hand side are 
the same graphs transformed according to the relation inverse 
to Eq. (5). 



space. We include in it the energy-momentum conserva- 
tion 5-functions. But we do not include there the rela- 
tivistic l/uj(ki) "on-shell" factors, which instead appear 
in the "matrix element" of the corresponding process, 
Uk- This will make the discussion of relativistic and non- 
relativistic cases uniform. The function F[n] is a sum of 
products of the type n^ 1 YYiLi n ki , where j S {1, . . . , m} 
with appropriate signs and permutations of indices for in- 
coming and outgoing particles. All dynamical aspects of 
turbulence follow from the scaling properties of the sys- 
tem [15]. Let u>k,nk and Uk have defined weights under 
a ^-rescaling of Fourier-space, 



.^k m ) = ^U(kt,. 



(7) 



The weight of the full collision integral under this re- 
parametrization is 



I, 



£d(ro-2)-a+/3+(m— l)7/ fe 



(8) 



It follows that the stationary turbulence with constant 
energy flux over momentum space is characterised by 
a power-law distribution function, n k ~ k~ s , where 
s = d + f3/(m — 1). The scaling properties also give the 
exponents of the self-similar distribution, Eq. (5). As- 
suming energy conservation in particles and with £ = t~ p 
we find q — Ap and p — l/((m — l)a — (3). For stationary 
turbulence we find that p should be (to — 1) times larger. 

For a A</> -theory in three spatial dimensions and four- 
particle interaction we have m = 4, (3 = —4a and a = 1. 
In this case s = 5/3 and p = 1/7. For three-particle 
interaction (the fourth particle belongs to the condensate 
in this case and the matrix element contains an additional 
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factor of (f> ) we find s — 3/2 and a smaller value for p 
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compared to the previous case. We can not distinguish 
between 5/3 and 3/2 for s in our numerical integrations, 
s rather fluctuates between these two numbers, while 1/7 
for p gives a fit to the data not as good as displayed in 
Fig. 3. However, during the integration time the energy 
in particles is neither conserving, nor there is a stationary 
source of energy. Namely, starting from the time at which 
the solution becomes self-similar, r <~ 3000, to the end of 
our integration, the energy influx from the zero mode to 
particles is about 20 %. Correcting for this energy influx 
we find q « 3.5p and p~ 1/6. This should be considered 
as satisfactory agreement given the simplifications which 
were made. 

Equilibration time and temperature. At late times the 
influence of the zero-mode should become negligible, but 
we still may expect the self-similar character of the evo- 
lution. Solution Eq. (5) with p = 1/7 should be valid 
in this case. This allows us to find the time needed to 
reach equilibrium. Indeed, the classical evolution will 
continue until the occupation numbers in the region of 
the peak in Fig. 3 will become of order one. At this time 
quantum effects become important and the distribution 
relaxes to thermal. Values of momenta were this happens 
are k max <~ A 1 / 4 ^^). On the other hand the initial dis- 
tribution is centred around fco ~ A 1 / 2 <po(''"o) and moves 
to ultraviolet according to Eq. (5) as cx k T p . It follows 
that the time to reach equilibrium is r ~ A~ 7 / 4 <~ 10 23 , 
where in the second equality we assumed the normaliza- 
tion to the inflationary model. For the reheating temper- 
ature we find, rotating back from the conformal reference 
frame, T R ~ k max /a(r) ~ A 2 po(t ) ~ 10- 26 M P1 ~ 100 
eV, where for the conformal scale factor we have used 
a(r) = t. 

Conclusions. Reheating after preheating appears to 
be a rather slow process. Although the "effective tem- 
perature" measured at low momentum modes during pre- 
heating may be high, in the model we have considered the 
resulting true temperature is parametrically the same as 
what could have been obtained in "naive" perturbation 
theory. Namely, equating the rate of scattering in ther- 
mal equilibrium to the Hubble expansion rate one ob- 
tains T ~ A 2 Mpi in this model. We anticipate this result 
should be applicable to more realistic models of inflation. 
Note that realistic models involve many fields and inter- 
actions and larger coupling constants will determine the 
true temperature. 
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